Ni coarsening in the three-phase soHd oxide fuel cell anode 
- a phase-field simulation study 



Hsun-Yi Chen% Hui-Chia Yu^ J. Scott Cronin'', James R. Wilson'', Scott A. Barnett'', 

Katsuyo Thornton'' 

"Department of Materials Science and Engineering, University of Michigan, 2300 Hayward St, Ann Arbor, 

MI 48109, USA 

''Department of Materials Science and Engineering, Northwestern University, 2220 Campus Drive, 

Evanston, IL 60201, USA 



Abstract 

Ni coarsening in Ni-yttria stabilized zirconia (YSZ) solid oxide fuel cell anodes is 
considered a major reason for anode degradation. We present a predictive, quanta- 
tive modeling framework based on the phase-field approach to systematically examine 
coarsening kinetics in such anodes. The initial structures for simulations are experi- 
mentally acquired functional layers of anodes. Sample size effects and error analysis 
of contact angles are examined. Three phase boundary (TPB) lengths and Ni surface 
areas are quantatively identified on the basis of the active, dead-end, and isolated phase 
clusters throughout coarsening. Tortuosity evolution of the pores is also investigated. 
We find that phase clusters with larger characteristic length evolve slower than those 
with smaller length scales. As a result, coarsening has small positive effects on trans- 
port, and impacts less on the active Ni surface area than the total counter part. TPBs, 
however, are found to be sensitive to local morphological features and are only indi- 
rectly correlated to the evolution kinetics of the Ni phase. 

Keywords: solid oxide fuel cell, coarsening, phase-field model, three-phase 
boundary, nickel, yttria-stabilized zirconia 



1. Introduction 

1.1. Nickel coarsening in SOFC anodes 

Solid oxide fuel cells (SOFCs) are one of the most promising clean energy con- 
version devices for stationary applications because of their low pollutant emissions, 
high efficiency, and abiUty to operate using various hydrocarbon fuels. The need for 
precious-metal catalysts is eliminated in SOFCs because the reaction kinetics is en- 
hanced at their operating temperatures, which are between 500 and 1000°C iH. How- 
ever, high operating temperatures also lead to disadvantages such as slow startup, high 
fabrication costs, and rapid component degradation |01. For SOFCs to be commercially 
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viable for stationary applications, their lifetime must meet a minimum of ~50,000 
hours. This goal has not yet been achieved. Understanding SOFC degradation mecha- 
nisms is therefore crucial to improving their durability. 

The degradation mechanisms of SOFCs have been reviewed in a few articles ^ 
Ql- SOFC degradation is usually evaluated in terms of cell power or voltage de- 
crease or area-specific resistance (ASR) increase, which can be determined by the 
electrochemical impedance spectroscopy (EIS). These electrochemical methods al- 
low the cell degradation to be monitored in situ. However, even though a few stud- 
ies isl @] have attempted to address SOFC degradation, it is very difficult to unam- 
biguously characterize the impacts of individual mechanisms using these techniques. 
Many non-electrochemical methods have been used to identify and monitor the degra- 
dation and failure mechanisms of SOFCs in situ or pre-/post-operation (for a review, 
see Ref. 01), including X-ray tomography and X-ray diffraction. Among these degra- 
dation mechanisms, the microstructural change in a Ni-based cermet anode is one of 
the least understood, because experiments that can provide detailed time-dependent, 
three-dimensional (3D) structural information are difficult. 

The anode of SOFCs is usually made of a composite material of complex mor- 
phology that facilitates electrochemical reactions, which require simultaeous transport 
of fuels, ions, and electrons. The most commonly used anode material up to date is a 
porous cermet comprised of Ni and YSZ. Electrochemical reactions at the SOFC anode 
mainly occur at the vicinity of TPBs, where the pore, Ni, and YSZ phases are in con- 
tact 181]. Although the electrochemical reaction mechanisms are not very well under- 
stood, the length of TPBs is considered one of the most important geometrical parame- 
ters that dictates the resistance in a SOFC anode 101 ■ The anode is thus designed to have 
a complex microstructure and increase the TPB lengths and the simultaneous transport 
of three different species. However, this intricate anode microstructure is typically not 
stable. The agglomeration and coarsening of nickel particles have been considered the 
major mechanisms responsible for microstructural change in SOFC anodes 113]. Ni 
coarsening in the SOFC anode is a capillarity-driven phenomenon. Regions with high 
curvatures have higher chemical potentials than those with lower curvatures in accor- 
dance with the Gibbs-Thomson effect. Materials will therefore be transported from 
higher- to lower-curvature regions when mobility is sufficiently large for the time scale 
of interest. This phenomenon increases the resistance in SOFCs and results in cell 
degradation. 

Long-term coarsening experiments have been conducted for thousands of hours to 
study Ni coarsening in Ni-YSZ anodes. Simwonis et al. ifioll measured a 33% decrease 
in electrical conductivity after 4000-hour exposure of a Ni-YSZ anode to a H2 environ- 
ment at 1000°C. They also found a 26% average Ni particle-size increase via analysis 
of micrographs of cross-sections. Thyden et al. jsllTll] performed an aging experiment 
of a SOFC over 17,500 hours. The cell was operated at 850°C with an initial current 
density of 1 A/cm-. Optical microscopy, field emission-scanning electron microscopy 
(FE-SEM), SEM-charge contrast (SEM-CC), focused ion beam (FIB)-SEM, and EIS 
measurements were utilized to analyze the microstructural evolution in the Ni-YSZ an- 
ode. Their results suggested that an increase in the H2O concentration can promote Ni 
particle coarsening and lead to conductivity loss within the Ni-YSZ cermet. Tanasini 
et al. I12II also conducted coarsening experiments for a single SOFC operated at 850°C 
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with humidified H2. Ahhough cell degradation is often attributed to cathodic processes 
(e.g., Cr poisoning and cathode-electrolyte interface formation reaction), they reported 
that the major cell performance reduction stems from the anode degradation due to 
coarsening. They also found that the cell potential drop and Ni particle size increase 
reach a plateau after ~ 1000 hours of operation. Despite these experimental efforts, a 
quantitative correlation between microstructure and coarsening still awaits explication. 

Even though the performance of SOFCs can be largely affected by microstructural 
changes, only a few models have been proposed to study the effects of coarsening on the 



performance of electrodes II13I - I15I1 . Due to the inability of acquiring 3D microstructural 
information for Ni-YSZ anodes, these models were mostly based on empirically fitted 
parameters or simplified microstructures. The use of oversimplified microstructures 
and empirical parameters without extensive validations can be problematic because it 
can lead to incorrect conclusions. 

Recently, we have demonstrated that the microstructure of SOFC electrodes can 



be three-dimensionally reconstructed using dual-beam FIB-SEM III6II . The data can 
be acquired to produce a wealth of information, including the tortuosity of individ- 
ual phases and the TPB density. While FIB-SEM allows the reconstruction of anode 
microstructure in three dimensions, this technique damages the materials. Thus, no 
evolution information can be obtained for the specimen after the procedure. Model- 
ing offers the advantage of allowing systematic analyses of coarsening effects on the 
performance of SOFC anodes. 

The phenomenon of anode coarsening can be described as a free-boundary problem 
in the sharp-interface modeling framework; however, explicit tracking of the evolving 
phase boundaries is highly impracticable in three dimensions. For a 3D anode with 
a complex microstructure, it is therefore advantageous to use a diffuse interface ap- 
proach, such as a phase-field modeling, to model the microstructural evolution. There- 
fore, we have developed a diffuse-interface modeling framework, based on the phase- 
field model and the smoothed boundary method (SBM), in conjunction with FIB-SEM 
experiments to quantitatively investigate Ni coarsening. 

1.2. The phase field model 

The fundamental basis of a phase-field approach is to define field variables, or 
order parameters (OPs), that distinguish the phases in a multi-phase system. OPs have 
a constant value in each bulk phase, while interfaces are represented by finite regions 
in which the OPs smoothly vary from one bulk value to another Since the boundary 
information is embedded in OPs, explicit tracking of moving boundaries is no longer 
necessary. This leads to a computational advantage in modeling multi-phase systems 
with multiple interfaces, especially in three dimensions. The phase-field model, one of 
the phase-field approachs, can also be considered as a type of diffuse interface model, 
which describes interfaces using a finite thickness. 

The inclusion of OPs into the free energy density is generally attributed to Lan- 
dau and Ginzburg for their work related to the modeling of superconductivity in the 
1950s UtIi . Later, in order to describe the interfacial energy of an inhomogeneous sys- 



tem, Cahn and Hilliard HI8II proposed to describe the free energy of a system by its OPs 
and their spatial derivatives. The concept of describing the evolution of an interface 
between two phases differing in composition with a Ginzburg-Landau-type functional 
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was introduced by Langer ifioll . Allen and Cahn ||20[l developed the theory for the mo- 
tion of a coherent anti-phase boundary. Formal asymptotic analyses have been used 
to show that a variety of phase-field models (PFMs) recover the corresponding sharp 
interface models when the width of the interfaces approaches zero fl2lll22ll . 

PFMs have been successfully utilized to descri be p hase transition in two-phase 
systems (see, e.g., 12111 ). dendrite growth (see, e.g., 12311 ), and Ostwald ripening (see, 
e.g., II24II ). However, in three-phase systems, a PFM with a generic free energy func- 
tional may introduce an artificial third-phase contribution at a two-phase boundary, if 



no special treatment is applied 125L 12611 . This effect can compromise the study of cer- 
met anode coarsening since the emergence of a phantom phase at a two-phase boundary 
can contribute to extra TPBs. This foreign-phase creation will lead to an erroneous es- 
timation of the electrochemically active sites. Nestler et al. l25ll proposed two types 
of remedies for this problem in a multi-component liquid-solid system. In each of 
these approaches, the potentials penalize equal contributions from each of the OPs to 
reduce the third phase appearance. However, the modified potentials may not eliminate 
the PFM simulation artifact for a simple three-phase conserved-field system undergo- 
ing coarsening. To address this issue, Folch et al. u.m developed a specific minimal 
model for a three-phase system that ensures that there is no third phase invasion at a 
pure two-phase boundary. This model can also accommodate systems with unequal 
surface tensions by adding tunable saddle-lifting terms. However, the computations 
become too expensive due to the steeper free-energy landscape needed for large-scale 
simulations. 

The PFM handles the interface implicitly because the interface information is em- 
bedded in OPs. Immiscibility and interfacial energy are naturally incorporated in 
the PFM via the bulk free energy and the gradient energy penalty over the difFuse- 
interfacial region. The interfacial energy ratios among diff'erent interfaces can be spec- 
ified in the PFM using proper parameterization of the free energy functional. There- 
fore, we developed a first phase-field model, named Model A, within the context of a 
multiphase PFM, to model the Ni coarsening in three-phase anodes i27ll . As discussed 
later. Model A allows a small mobility of YSZ, which significantly aff'ects the evolu- 
tion. Therefore, we developed an alternative model based on the smoothed boundary 
method (SBM), which is briefly described in the following section. 



1.3. The smoothed boundary method 

Complex geometries are abundant in naturally and man-made objects. To study 
physical processes or phenomena occurring within such objects, the numerical solution 
of partial differential equations (PDEs) with prescribed boundary conditions is neces- 
sary. A standard scheme requires triangulation of the complex shapes, followed by 
solving of the PDEs using the finite element method (FEM). However, automatic gen- 
eration of proper meshes for complex 3D domains is challenging. In addition, in many 
cases the complex geometries can evolve during a physical process, thus demanding a 
dynamic remeshing of the evolved domain. 

To address these difficulties, one can alternatively embed complex geometries within 
a larger, simpler domain (such as a cube). The PDEs can then be solved with regularly 
shaped meshes on the extended domain, provided that the original boundary condi- 
tions can be properly applied. Methods employing such a concept include composite 
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FEM 112811 . extended FEM 112911 . the immersed interface method yOfl, the immersed 
boundary method |31], and the cut-cell method js^]- 

The smoothed boundary method differs from the above mentioned methods in that 
it represents the complex geometries with a phase-field-like function, which veries 
smoothly across the domain boundary. Thus, the sharp boundaries of the complex 
geometry dO. are instead described by a thin interface of a finite width. This diffuse- 
domain approach was utilized to study the diffusion of chemoattractant inside a cell 
with a no-flux boundary condition (BC) at the cell surface j fssll . Spectral methods were 
later coupled with the SBM to model electrical wave propagation in cardiac tissues 
with a no-flux BC HSl. In Ref. ^36^, the SBM was extended to solve PDEs in 
complex geometries with Dirichlet, Neumann, and Robin BCs. An alternative but 
similar approach was independently developed, |37:|. 

In the SBM, an auxiliary variable, a domain parameter iff, is introduced to identify 
the domain of interest Q. in which the PDEs are solved. The domain parameter i// 
usually has a value of 1 within Q, exterior to Q, and < i/' < 1 in dO.. Since the 
complex geometry is embedded in a regularly shaped, expanded domain, the exterior 
of the domain of interest is included only to facilitate computation; any numerical 
solutions obtained in the external region are devoid of meaningful information. The 
initial construction of the domain parameter can be achieved by solving a phase- 
field equation or by transforming the distance function of the domain structure with the 
hyperbolic tangent. 

It is generally believed that the YSZ phase has a very low mobility in an operating 
SOFC anode; that is, YSZ serves as the supporting structure in which Ni coarsens. 
We therefore propose a model, named Model B, utilizing the SBM for Ni coarsening 
simulations. In Model B, we assume that YSZ is stationary. It is therefore treated as 
the geometry within which the Ni and pore phases evolve. We also assume that the 
triple junctions in the Ni-YSZ anode possess the contact angles deduced from Young's 
equation for a locally flat surface. The dynamics of this system can thus be modeled by 
a single OP Cahn-Hilliard equation with two complementary BCs implemented with 
the SBM: the contact-angle BC at triple junctions and the no-flux BC at YSZ interfaces. 



2. Methods 

In this section, we elaborate the model framework that we previously published in a 



short communication for Ni coarsening ll27ll . The non-dimensionalization, asymptotic 



analysis, and numerical methods are discussed in detail. Moreover, the methodology 
we used to identify the TPBs and distinguish the percolated or isolated phases is dis- 
cussed. An error analysis is performed to examine how the ratios between the charac- 
teristic length and the interfacial thickness affect the contact angles at triple junctions 
in our model. 

2.1. Model Formulation 

We have proposed two PFMs for coarsening simulations in three-phase SOFC cer- 
met anodes. Our models are based on a free energy funtional that is computationally 
inexpensive and can circumvent the phantom phase issue associated with some other 
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solutions. The approximations we made in both of our models are: (1) the volume 
(mass) of each phase in the Ni-YSZ anode is conserved so that the Cahn-Hilliard dy- 
namics is applicable; (2) material properties required for modeling are estimated at 
1000°C; (3) Ni surface free energy is assumed to be isotropic; (4) Ni surface diffusivity 
is estimated by the mean value of multiple crystal orientations; (5) surface diffusion 
dominates. 

In our PFMs, the Cahn-Hilliard evolution equations can be considered as a type 
of diffusion equation, where the material transport is driven by the chemical poten- 
tial. This potential can be formulated by the design of a free energy functional that 
is based on the effective driving forces for material transport. At ordinary operating 
temperatures, the three phases comprise the Ni-YSZ anode can be assumed immisci- 
ble; therefore, lattice misfit as well as elastic energy contribution toward the diffusion 
are negligible. In addition, external forces such as gravity are not considered. We thus 
formulate our free energy functional as the simplest Ginzburg-Landau type functional 
F - J fdV, where the free-energy density / simply depends on the OPs and their 

V 

derivatives. 

The two PFMs we have proposed differ in their treatment of the YSZ phase. If we 
consider that the YSZ phase can transport with its mobility, the Ni-YSZ anode is then a 
system of three mobile phases (Ni, YSZ and Pore). This system requires two evolution 
equations to resolve the kinetics. We call this model Model A. In Model B, we consider 
the YSZ phase to be completely immobile, which is justified because its mobility is 
orders of magnitude smaller than that of the Ni phase. Therefore, only two phases, 
Ni and pore, are allowed to evolve in between of the YSZ matrix. To incorporate the 
YSZ phase as the internal boundary in the computational domain, we have developed 
a model utilizing the smoothed boundary method. 

2.7.7. Model A 

In Model A, the OPs are vectorized as (p - {(f)i,(p2,4>i) for generalization. Each 
vector component represents the volume fraction of the corresponding phase so that an 
additional constraint 0i + 02 + 03 = 1 applies; thus, only two evolution equations are 
needed. The governing equation sets can be written as 

^ = V-M(0)V— , 
ot 0<p\ 

(1) 

^ = V-M(0)V— , 

ot 002 

where M( ) is the surface mobility. 

A standard bulk free energy for a three-phase system should contain three local 
minima, each of which represents a bulk-phase value. We select the three local minima 
as Ofl = (1,0,0),<E)/, = (0, l,0),<l)c = (0,0,1). However, in the three-phase Cahn- 
Hilliard dynamics, a foreign phase can be introduced at a two-phase interface. Two 
remedies are that one can use the procedure developed in ll26ll to generate a specific 



free energy, or one can use an interpolation function similar to that described in [38 1 



(the latter method may only reduce the amount of the third phase contribution). Since 
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the three phases in the Ni-YSZ anodes are immiscible, the generation of a foreign phase 
at a two-phase boundary is a problematic artifact that leads to extra TPB sites that do 
not physically exist. 

Our approach to resolve this issue is to select a free energy functional that can 
limit a foreign-phase appearance and can also incorporate an unequal surface tension 
at a reasonable computational expense. In the Cahn-Hilliard dynamics, the excess free 
energy at the interfaces, i.e., the interfacial energy, is determined by gradient energy co- 
efficients and the bulk free energy at intermediate values of OPs. According to Young's 
relation, the contact angles among the different phases are determined by the interfacial 
energies of the intersecting interfaces. For instance, if phases 1, 2, and 3 are in contact 
at a triple junction, the contact angle formed in phase 1 is related to the interfacial en- 
ergy of the interface 2-3 relative to that of the interfaces 1-3 and 1-2. It is, however, 
easier to understand the interfacial energy as a factor reflecting the affinity of one phase 
to another. A lower interfacial energy suggests a stronger bonding of the two phases 
that are in contact. This tendency appears to have profound effects on the long-term 



microstructural evolution as demonstrated in our simulations II27I1 . 

Applying the commonly used gradient and bulk free energy densities for multiple- 
OP systems, our free energy functional can be written as 



r ( ^ 3 



3 

(2) 



\i=\ i.j=l,i>j 

The specific interfacial energy is calculated from the equilibrium solution at a pla- 
nar interface. To do so, the free energy functional must be minimized subject to a con- 
straint that the sum of the order parameters equals unity. This is carried out by taking 
the variational derivative of the functional using the method of Lagrange multipliers. 
The resulting equation is given by 



6F\ 



- 



_6F SF 



where the variational derivatives at the right hand side of the equation are taken as if 
all <^,'s were independent ll26ll . If we consider an interface between phase / and phase 
i without the existence of the third phase, the interfacial free energy can be calculated 
from the integration of the free energy density over the interface as 



jij = 2aij ^{ki + kj) j| 0/( 1 - <^/)#/ = Y ^j(k~7k~). (4) 
Similarly, the interfacial width can also be found as 



an 

6ij = \ (5) 



After multiplying Eq.|4]by Eq.|5]and rearranging the result, we find 



35ijyij = a^. (6) 
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By choosing Sn - S23 - ^13, the relationships between the gradient energy coefficients 
and the interfacial energy in this three-phase system is simpHfied as 

2 2 2 ' ^ ' 

a^i '^n 

The values of the surface free energies can be set to those found in the literature or be 
left as a parameter if data are missing or uncertain. The interfacial width is a computa- 
tional parameter that must be chosen to ensure thin-interface limit for accuracy while 
providing a numerically well-resolved interface. Once the interfacial width and in- 
terfacial energies are set, the gradient energy coefficients Oij can be calculated from 
Eq.|6l and the bulk energy coefficients A:, can be calculated from Eq.|5] 

Surface diffusion is generally considered to be the main mechanism of transport for 
the capillarity-driven microstructural evolution in these SOFC anodes ifTsll . We have 
confirmed the validity of this assumption based on an analysis of order of magnitudes. 
Since the OPs with integer values represent the bulk phases in the PFMs, the surface 
mobility function in a one-OP, two-phase case is commonly formulated as M((p) - 
0^(1 - (p)^\ such a formulation guarantees a non-zero mobility at two-phase boundaries 
only. 

In Model A, the surface mobility function is given by: 

3 

, 02, 03) = Yu ^'j nc,c„ (<^<) ncQ. (<PM<Pj( 1 - 1 - 0;)), (8) 

'.j=iJ>j 

where we introduce a boxcar function nc/C,,(0() to avoid excess mobility resulting 
from the appearance of a small amount (less than 3% in fraction) of a foreign phase 
at a two-phase boundary. We choose C; - 0.05 as the lower cutoff OP value in the 
mobility function, and C/, = 1 - C; as the upper cutoff value and resulting from the 
complementary value of the lower cutoff. All three mobility prefactors have positive 
values, and thus, materials flow from high to low chemical potentials. 

2.7.2. Model B 

In this model, we assume that YSZ is stationary and that the Ni phase evolves 
by diffusion along the Ni-pore interfaces. The dynamics of this system can thus be 
described by a single-OP Cahn-Hilliard equation with one contact-angle BC at triple 
junctions, and one no-flux BC at YSZ interfaces, where the OP distinguishes the Ni 
and pore phases. 

The treatment of a contact-angle BC in a non-conserved (Allen-Cahn) PFM was 



proposed by Warren et al. 113911 to model heterogeneous nucleation. Following [SS 
we developed Model B using the SBM for a Ni-YSZ anode coarsening with conserved 
dynamics. This SBM approach is an entirely diffuse-interface treatment that implicitly 
handles complex geometries (while Warren et al. applied delta function ^^\). In the 
following derivation, we demonstrate that both no-flux and contact-angle BCs can be 
coupled with a Cahn-Hilliard equation to give a single evolution equation using the 
SBM framework. 
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For one-OP Cahn-Hilliard dynamics within a SBM framework, the evolution equa- 
tion can be written as 

§ =V-M(0,(A)V//. (9) 
ot 

We consider the free energy functional as being: 

5" = j^flfV[y|V^l' + /(0)], (10) 

where /(0) is a generic double-well function. The chemical potential is by definition 
the variational derivative of free energy T with respect to the order parameter ^, i.e., 
jj - Sf/Sip - df ld(p - e^V^0. We introduce in the SBM a domain parameter i/r to 
incorporate general BCs on internal boundaries. In this case, the domain parameter ^ 
distinguishes regions having YSZ (i/r = 0) and other phases {ijj - 1), as well as the YSZ 
interfaces (0 < i/' < 1). 

Multiplying Eq. (|9ll by (/r and letting J - MVfi the formulation becomes 

-^^^W-J^V-(ifrJ)-J-Wi(r. (11) 
ot 

The aforementioned no-flux BC, Vju ■ Vi/* = 0, should be applied to the internal bound- 
aries (YSZ interfaces) to ensure mass conservation. This no-flux BC eliminates the last 
term of Eq. (fTTI) and the equation becomes 

^ V ■ (^7l = V ■ (iJfMVfi). (12) 

Ot 

Like an ordinary OP in a phase-field approach, varies continuously across the 
interface; thus, the unit interface normal ft of the YSZ interface can be described as a 
function of the gradient of i//, i.e., it - Vi/// By assuming that the effects of the 
YSZ on the Ni phase are only of short range (immiscible) and <p = I represents the 
bulk Ni phase, the contact angle 9 at triple junctions can be formulated as 

^ V0 ViA V(* 

n ■ = ■ = - cost/. (13) 

m m m 

The negative sign comes from the convention that Vifr points into YSZ and V0 points 
out of Ni. 

The mechanical equilibrium at the triple junction corresponds to an extremum of 
the free energy, i.e., d'F - 0. We can use the planar solution of the thermodynamic 
equilibrium condition within the interface to find a useful equality, |V0| = ^jlfle, 
which can be substituted into Eq. ( fT3b to derive the SBM contact-angle BC as 

Vi/r ■ V0 = - I Vi/'l COS 9-!—- . (14) 
e 

This contact-angle BC results in energy only near triple junctions, rather than in the 
bulk volume, to achieve the force balance at the triple junction that is dictated by 
Young's equation for a flat surface. The final evolution equation derived from the SBM 
with no-flux and contact-angle BCs is: 
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dt 



V ■ {(AMV[/^ - -(V ■ + "^'J COS 0] 



, , , (15) 

Again, the gradient energy coefficient e, the contact angle 6, and the bulk energy co- 
efficients in /(0) are determined according to the interfacial energies of the Ni-YSZ 
cermet and the selected interfacial widths. 

The mobility function in Model B is formulated as 

M(cP, lA) = MNi-Por. nQC„(0)(<^'(l - <P^))8W, (16) 

where g{ip) = if/^{lOifP'- 15i/' + 6) is introduced to control the mobility at and near triple 
junctions. This one-sided interpolation function g((/') transitions smoothly from 1 to the 
order of 0.01 as the domain parameter varies from 1 to 0.5. In other words, this choice 
of the mobility function ensures an immobile YSZ phase by limiting the mobilities at 
YSZ interfaces. Subsequently, the mobility near a triple junction decreases from the 
Ni-pore value to a value that is 10"^ smaller as i/r varies from 1 to about 0.1. 

2.2. Nondimensionalization and Asymptotic Analysis 

Having the appropriate values of the mobility prefactors is essential in simulating 
coarsening kinetics. In order to quantitatively correlate simulation results with physical 
phenomena, asymptotic analyses are required in PFMs because of their diffuse inter- 
face nature. Using an asymptotic analysis, we determined the relationship between 
diffiisivity and mobility and the characteristic simulation time scale. Unlike some of 
the previous studies, which determined model parameters by fitting coarsening exper- 
imental results, our model is a predictive model that is free of fitting parameters when 
the material properties, such as surface diffusivities, are accurately known. 

In the Ni-YSZ anode of SOFCs, coarsening proceeds mostly via surface diffusion. 
The anisotropic effect of crystal facets on surface diffusion is lumped into one ensemble 
diffusivity in our models. In this case, the normal velocity V„ of the interface T incurred 
by surface diffusion, in the dimensional sharp-interface form, can be represented as 

kbTNv kbTNv os- 

, where is the surface energy, is the surface diffusion coefficient, 5,, is the inter- 
facial thickness, A^,, is the atomic number density per volume, is the local curvature, 

is the surface Laplacian and d/ds is the gradient operator along the interface 1I40I1 . 

To link the corresponding diffuse-interface model to the aforementioned sharp- 
interface model, a pure two-phase boundary in this three-phase system is considered. 
For both Model A and Model B, the dimensional Cahn-Hilliard evolution equation at a 
two-phase boundary can be re-arranged as 

^ = V ■ (M(0)V/.), (18) 
ot 

;.= ^-.V0, (19) 
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, where (p is the OP distinguishing two phases, /(0) = W4>^{1 - (f>)^/4- is the energy 
density, and M(0) = 6Ms0^(l - <p)^ is the surface mobihty. Eq.[T8]and Eg. [T9l possess 
the physical steady-state solution for a planar interface, i.e., <p{x) = [1 - tanh(x/25)]/2, 
with the interface thickness 6 = e V(2/W), and the interfacial energy y - e V(W72). 



Following the derivation in 114 lH . a new set of variables for a non-dimensionalization 
of the governing equations is utilized. We assume that r = L*/D, wherein L is the 
characteristic length scale of the sample and t is the characteristic time scale for surface 
diffusion, and consider that other nondimensional quantities (denoted by overbars) are 
defined by 

36M,yd -6 X _ t 

D L L T 

m) = cp\\ - 0)2, /(0) = - (20) 

Eqs. [18] and [19] can be written using the non-dimensional parameters to derive the 
nondimensional evolution equation as 

^ = lvM(0)V/z (21) 
= — V 0, (22) 

00 



The asymptotic analysis derivation is similar to 114211 and is detailed in the appendix. 
Comparing the dimension-restored equation derived from the asymptotic analysis with 
the corresponding sharp interface equation, Eq. [17] we find 

1 = M,8y = 1^_L_I, (23) 

, which indicates that = DsSslkBTN^.S, provided the choice of y = y,; that is, we 
find that the surface mobility is proportional to surface diffusivity. Choosing - 36, 
the model variable D is connected to the surface diffusivity Ds as follows: 

D = (24) 

The time scale that links the simulation time to the physical time is acquired from 
Eq.Hias 

ysDsOs 

2.3. Numerical methods 

One of the challenges in modeling Ni-YSZ anode coarsening is the fact that solving 
Cahn-Hilliard equations with an explicit time iteration scheme is too expensive, espe- 
cially for a large scale simulations in three dimensions. Thus, we solve the nondimen- 
sional evolution equations, Eqs. [21] and [22] with the algorithm based on splitting the 
fourth-order Cahn-Hilliard equation into two second-order equations and soloving for 
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the OPs and chemical potential simultaneously iRT]. We use the central-differencing 
method for the spatial discretization and the Crank-Nicholson algorithm for the time 
discretization. This semi-implicit scheme significantly reduces the stiffness of the nu- 
merical integration, which allows much larger time step size. 

To solve the nonlinear finite-difference equations, Newton's method is utilized 
for the nonlinear terms. A pointwise Gauss-Sidel relaxation scheme and a red-black 
checkerboard iteration scheme are used together to accerelate the convergence rate as 
well to facilitate the parallelization. In most cases, we find this solver is over 100 times 
faster in comparison to the explicit scheme. However, attention must be paid to the 
choice of the time stepping size. We find in some cases, when an overly large time step 
is used in our solver, the simulation results are incorrect even if the numerical scheme 
is stable. 

When solving the coupled governing equations in Model A, the solver is paral- 
lelized with Message Passing Interface (MPl) library to take advantage of the multiple 
processors. In our MPl code, the domain is decomposed equally in size in each axial 
direction, if possible, to achieve load balance. For instance, if 64 CPUs are allocated 
to the computation, the entire domain is decomposed into 4 by 4 by 4 sub-domains; 
that is, each sub-domain has a domain size of 1/64 of the original domain in volume or, 
more explicitly, 1/4 of the length of the original domain in each axis (assuming there 
are no residuals). 

The solver for Model B is parallelized with both the MPl and openMP libraries. 
Since only one Cahn-Hilliard equation is solved, it is much more numerically stable 
and efficient compare to the Model A solver 

2.4. Error analysis: contact angles and the interfacial width 

The phase-field model is known to smooth out microstructures with length scales 
below the diffuse interface thickness. This artificial smoothing process introduces er- 
rors during the early stages of a simulation, but has no negative effects on the analysis of 
long-term coarsening kinetics. In contrast, microstructural evolution kinetics can only 
be accurately resolved in PFMs when the length scale of the microstructure is larger 
than certain multiples of the interfacial thickness. In other words, there is a critical ratio 
between the microstructural length scale and the interfacial thickness that is required 
to obtain a sufficient agreement with the sharp-interface limit. For example, if a system 
contains particles with typical radii a few times smaller than the interface thickness, the 
simulations of its evolution will incur a large error The commonly recognized critical 
ratio is of about 10. However, in large scale 3D simulations, to resolve a complicated 
system with microstructural features of different length scales, achieving this critical 
value for all features is impracticable because it would require a very high resolution. 
In practice, we use a smaller ratio, especially when there are smaller features among a 
range of feature sizes within the microstructure, and quantify the errors introduced by 
the selected value. 

Because the coarsening kinetics of Ni- YSZ anodes has been found to depend strongly 
upon the Ni-YSZ contact angle, an error analysis is performed via the investigation of 
the contact angle of Ni on the YSZ phase in two dimensions (2D). The two-dimensional 
domain is initialized with a bottom that is half occupied by YSZ and a top half that is 
equally divided into the Ni and pore phases with a 90° contact angle. The remainder 
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of the domain is filled with the pore phase. The domain size is designed to be large 
enough so that no-flux BCs have negligible effects on the contact angles at the triple 
junction. The system is evolved with our Model B to its steady state. The final contact 
angle 6c is calculated from the average value of the dot product of the normal, 

■ 7;^ = - cos (26) 
IV-AI \^<P\ 

over the region where the domain parameter i/r and the order parameter 4> are both 
between 0. 1 and 0.9, indicating the TPB region. 

Two contact angles are studied: 120° and 93°. The 120° case is selected as a 
reference case because the cosine function is away from the extrema or inflection point 
at this contact angle. The 93° case corresponds to a physical contact angle of Ni on the 
YSZ phase that is based on our selected interfacial energies. 

As the characteristic length of the system, we choose the domain size, which is 
varied from 10 points to 100 points in each direction, while the interfacial width is held 
at 4 grid points. Therefore, the ratio of the domain size to the interfacial width varies 
from 2.5 to 10. For the 10 x 10 domain, the interfacial region (either the domain or 
order parameters are between 0.1 and 0.9) occupies about half of the domain. As the 
domain size increases, the fraction becomes very small. The contact angles normalized 
to the set value versus the ratios of the domain size to the interfacial width are plotted 
in Fig. [1] In the case of a 120° contact angle, we find that the angle deviate from the 
sharp-interface value by less than 0.2° when the ratio is larger than 5. Even at a ratio 
of 2.5, the deviation is below 1 degree or 1%. In the case of a 93° contact angle, the 
contact angle differs from the sharp-interface value by less than 0.5°, or 0.5%. 
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Figure 1 : Normalized angle plotted with the ratio of the domain size to the interfacial width. The error 
decreases with an increasing ratio of the domain size to the interfacial width. 
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The interfacial width utilized in coarsening simulations is approximately 0.13 j-im. 
Based on our error analysis results, the contact angles of the Ni particles in the active 
and dead-end categories with a length scale larger than 0.32 /im can fairly accurately 
be modeled by using our Model B. In the isolated network, a length scale larger than 
0.64 jum is needed for a fair accuracy. For the (4 //m)^ domain, only 3% in volume 
of the Ni phase belongs to clusters outside the accurate range. In addition, based on 
Young's equation for a flat surface, a 0.5° difference in contact angle of the Ni phase 
corresponds to a 1.1% difference in Ni-YSZ interfacial energy. In other words, this 
model error has very limited effects on the coarsening kinetics of the Ni-YSZ anode. 



3. Parameterization of the models 

The aforementioned asymptotic analysis provides the mathematical ground for cor- 
relating simulations on the coarsening phenomenon with experimental results. How- 
ever, to quantitatively model the physical process, the material-specific parameters in 
the governing equations must be specified based on the material properties. Most of 
these parameters can be found in the literature. 

3.1. Mobility prefactors 

As demonstrated in our asymptotic analysis, the mobility prefactor correspond- 
ing to the Ni-pore interface should be proportional to the surface self-diffusivity of 
Ni, which is usually anisotropic or depends on the crystallographic orientation. The 
Ni diff'usivities based on the field ion microscope (FIM) measurements range from 



10"' ' m-s-' in (110) to 10"'' m^s"' in (331) at 1273 K |43M45|]: these values are ex- 
trapolated at high temperatures and could be inaccurate, because FIM can be conducted 
only at low temperature regions (T < Q.2T,„, T,,, is the melting temperature). On the 
other hand, the surface smoothing method (SSm) measures mass-transfer diff'usion at 
high temperatures (T > 0.77,,,), which is generally averaged over crystallographic ori- 
entations f43\. The latter method provides a more appropriate value for the diffusivity 
(~ 10"'" m^s"') during Ni coarsening at the modeled temperature (1000°C). In addi- 
tion, the anisotropy of the Ni surface diffusivity was found to be relatively small in high 
temperature regions (0.82 T„, < T < T,„) in |46|]. We thus consider an ensemble value 
~ 10""^ m^s"' for the Ni surface diffusivity in our model. 

Cation diffusion in oxides is related to various parameters, such as valence, atomic 



radius, impurity, and oxygen activity. In Ref. 114711 . it is reported that the bulk diffusivity 
of yttrium in YSZ (containing 10 to 32 mol% Y2O3) is slightly larger than that of 
zirconium in YSZ, and that the difference is less than an order of magnitude. The 
surface diffusion of Zr in YSZ is calculated from the measurements of the surface area 



reduction in powder compacts during sintering in Ref. Il48ll . and the diffusivity at 1000 
°C is found to be ~ 10"'^ m^s"'. We thus consider a surface diffusivity of YSZ of 
10"'^ m^s"' in our simulations. 

The diffusion mechanisms at metal-ceramic interfaces are less understood than 
those on metal or ceramic surfaces. The cohesiveness of the Ni-YSZ interface, which 
depends on the process used to fabricate the porous cermets, determines the interfacial 
structure and affects the effective interfacial diffusivity. A series of diffusion bonding 
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experiments indicated that the metal-ceramic interface does not act as an efficient va- 
cancy sink or mass transport path [49]. These findings suggest that the difFusivity at 
the Ni-YSZ interface should be much smaller than that of the Ni surface and the YSZ 
surface. In addition, the exact value of the diffusivity at the Ni-YSZ interface is not 
as important as that of Ni because the redistribution of Ni is observed to be the dom- 
inant morphological change during anode coarsening. The difFusivity of the Ni-YSZ 
interface is considered to be approximately ~ 10"^° m-s"'. 

Using these values, the mobility ratios among the three materials are set at Mnip : 
Myp : MNiY - 1 : 10"^ : 10"'°, where the subscript Ni represents nickel, Y the YSZ 
and P the pore phase. For non-dimensionalization, a mobility scale of 10 '° m-s ' is 
used. 



3.2. Bulk and gradient energy coefficients 

In a PFM, the balance between the bulk and the gradient energy terms in the free 
energy functional determines the thickness of the interface. The ratio of the interfacial 
thickness to the characteristic length of the system and the ratio of the domain size to 
the characteristic length are crucial for the validity of a PFM. It is known that an inter- 
face that is too thin can cause an unphysical pinning or halting of coarsening, whereas 
an interface that is too thick can lead to unphysical dissolution of particles jSOj]. There- 
fore, our task is to select bulk and gradient energy coefficients that result in an appro- 
priate interfacial energy, while keeping the interfacial thickness sufficiently large so 
that the computation is feasible. 

In a multi-phase system, using an unequal interfacial thickness among different in- 
terfaces without a special treatment can cause erroneous evolution kinetics in a PFM 
simulation because the interfacial thickness changes the ratio of the simulation time 
scale to the physical time scale. Therefore, a single value for the interfacial thickness 
of all interfaces is set. For example, we set 6ij = A.y, which gives 4 points in the inter- 
facial region to avoid the aforementioned pinning while optimizing the computational 
efficiency (A.y depends on the selected length scale in a non-dimensionalization). The 
interfacial energies are obtained from existing experimental or computational studies. 
These choices are then used in Eqs.|5]and|6]to determine the bulk and gradient energy 
coefficients. 

The three interfaces that exist in the Ni-YSZ anode are the Ni surface (Ni-pore 
interface), the YSZ surface (YSZ-pore interface), and the Ni-YSZ interface. The Ni 
surface free energy has been widely investigated and the values reported are in goo d 
agreement. At 1 000°C, the surface energy of ~ 1.9 Jm"- for Ni has been measured 15 ill . 
A 8 YSZ surface (8 mol % Y2O3) has been reported in Ref. fs^]. It was studied at 
1300~1600°C using a multiphase equilibrium technique. The YSZ surface free energy 
yYSz was found to decrease linearly from 1 .26 to 1 . 1 3 Jm - in the range of temperatures 
studied. By extrapolating this data, one obtains yyp ~ 1 -4 Jm"^ at 1000°C. This value is 
used to parameterize our model. In Ref. 15311 . ab initio calculations have been reported 
with yyp ranging from 1.04 to 1.75 Jm at T = 0, depending on the crystallographic 
orientation. 

The free energy of a heterogeneous Ni-YSZ interface depends on the interfacial 
structure, which is dependant on the fabrication process. Nikolopoulos et al. ex- 
perimentally measured the non-reactive contact angle between molten Ni and 8 YSZ 
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at 1500°C and found a value of 117°, which suggests an interfacial energy of 1.95 
Jm"^ I54II . Although this wetting experiment was not conducted at normal SOFC op- 
erating temperatures, the result implies a poor wettability at the Ni-YSZ interface. 

Several ab initio calculations on bond formation at Ni-YSZ interfaces have been re- 
ported ISJ . A strong bonding between Ni-Zr and Ni-O was found at the Ni( 1 00) /ZrOaC 1 00) 
polar interfaces, and an interfacial tension of cr(ioo) - 1 .04 Jm"- was reported. Another 
local minimum appeared at the Zr02/Ni(ll 1) interface where cr(iii) - 1.80 Jm"^ ll56ll . 
However, Ni-YSZ interfaces formed subsequent to sintering are multifaceted. Despite 
the lack of understanding of these interfaces, we know the range of reasonable interfa- 
cial energies at Ni-YSZ interfaces. In our simulations, yNiY = 1-50 Jm"^ is used. 

In view of the above discussion, we assume that the gradient energy coefficient 
ratios are yNiP : Q^Niv • '^yp = 1-9 : 1.5 : 1.4 based on a surface energy scale of 1 Jm"^, 
which leads to the bulk energy coefficient ratios of kj.^ : ky : kp = I : 0.5 : 0.9. 



4. 3D data analysis 

In three dimensions, a boundary at which three phases coincide has only one de- 
gree of freedom in space, i.e., the TPB is a line in the Ni-YSZ anode where the Ni, 
YSZ, and pore phases are all in contact. However, to identify the TPBs based on the 
microstructure reconstructed by FIB-SEM, some data processing is necessary since the 
microstructural data are represented by a 3D matrix. In this matrix, each phase is repre- 
sented by regions in which voxel values are equal to a pre-determined constant. There 
are several methods to determine the TPB length in 3D volumetric data. For example. 



in Ref. 115 711 . the TPBs are identified as the edges where three different voxels, each of 
which belongs to a different phase, are in contact. With some geometric corrections, 
the TPB length can be acquired with fair accuracy, but is limited by the resolution. 

In our work, we adopt a thinning, or skeletonization, algorithm to determine the 
TPB regions, and count the voxels. In PFM simulations, the interfacial region is com- 
monly identified as the zone over which an OP varies from 0.1 to 0.9. Because the 
interfaces in PFMs span multiple grid points, one can identify the diffuse TPB regions 
as voxels in which all three OPs are between 0.1 and 0.9. For Model B, the threshold 
values are slightly different due to the additional pre- and post- treatment of the data 
(Ref. [[ssf]). In order to recover the one-dimensional nature of TPBs from this data, 
we use a thinning algorithm developed in Ref. ll59ll . which reduces the diffuse TPB 
regions to the corresponding skeleton of chains of voxels. The key feature of this algo- 
rithm is that it preserves the topological characteristics of the original image/voxelated 
data (i.e., it does not allow the pinching or connection of regions). 

The final TPB length is calculated via a multiplication of the physical grid size and 
the total number of TPB voxels after skeletonization. The used procedures may either 
overestimate or underestimate the TPB lengths depending on the angle at which the 
TPB lies within the grid; for an isotropic distribution of lines, the over/underestimation 
has been determined to be less than 13%. Because we have nearly isotropic distribu- 
tion (i.e., anisotropy that only results from statistical variations), the error should be 
relatively consistent throughout the coarsening process and therefore not impact our 
investigation of the evolution of the TPB length. 
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Not all of the TPBs are active; active TPBs must be simultaneously in contact with 
the three active phases that facilitate simultaneous transport of fuel, electrons, and oxy- 
gen ions. Since the electrochemical reactions occur only at or near these active TPBs, it 
is these active TPBs that contributes to the anode performance. Before identifying the 
active TPBs, we must first identify the active phases. Physically, a Ni phase cluster is 
active only if it connects the TPB sites to the current collector; a YSZ cluster is active 
only if it connects the TPBs to the electrolyte; and a pore cluster is active only if it 
connects the TPBs to the gas flow channels. However, the determination of such con- 
nectivities is feasible only if an anode structure that spans from the electrolyte to the 
current collector and to the gas channel is available. Our reconstructed microstructure, 
in contrast, constitutes only a portion of the entire span. 

As a compromise, we follow the procedures described in Ref. 15711 . which cate- 
gorize each bulk region into active, dead-end, or isolated clusters. Active clusters are 
clusters that are connected to at least two sides of the sample domain boundaries, while 
isolated clusters are those that are not connected to any of the sides and are therefore 
electrochemically inactive. Dead-end clusters are defined as those that are in contact 
with only one side of the sample domain. In Ref. 15711 . a pre-smoothing procedure is 
used before this calculation to avoid artifacts resulting from image processing. Because 
PFM simulations naturally smooth microstructures over the lengh scale of the interfa- 
cial thickness, the pre-smoothing procedure is not necessary in our analysis. Note that 
in this procedure, a voxel is assumed to belong to a cluster if at least one of the 6 nearest 
neighbors is in the cluster The active TPBs are then identified as those simultaneously 
in contact with the active networks of all three different phases. The TPB is marked as 
inactive as long as one of the phases adjacent to a TPB voxel is isolated. The remainder 
of the TPBs are considered unknown in terms of activity. 



5. Results and Discussion 

5.7. Phase connectivity and TPB activity 

TPB length has been recognized as one of the most important geometric parameter 
in a three-phase SOFC anode that influences the electrochemical performance. How- 
ever, how the coarsening of microstructures affects the TPB activity has not been fully 
explored. As previously mentioned, TPBs are active only if they are simultaneously 
connected to conducting ionic, electronic, and gaseous transport pathways (phases). 
Therefore, the connectivity evolution of a phase induced by coarsening can alter the 
activity of TPBs and the effective conductivity of that phase. In addition, the evolu- 
tion of the anode microstructure results in changes in the amounts and distributions of 
TPBs. In turn, coarsening has significant impacts on the overall performance of the 
electrode. By simulating coarsening with our models, a series of microstructures are 
acquired at various stage of evolution. Using the methods described in section |4] we 
can analyze the evolution of the active TPB length based on the evolving TPB distri- 
bution and phase connectivity during coarsening. 

The initial microstructure of our simulation is based on an FIB-SEM reconstructed 



SOFC anode of dimensions 9.73 ^im x 8.35 jjm x 11.2 //m 111611 . We selected a 4.0 



//m X 4.0 /^m x 4.0 /^m portion of the specimen (Fig. |2]i, resolved by a domain of a 



17 




Figure 2: Initial Ni-YSZ anode microstmcture for a set of larger simulations. The dimensions of the mi- 
crostructure are 4 /im X 4 fim X 4 pm. The Ni, pore, and YSZ phases are represented in green, blue, and 
semi-transparent, respectively (with volume fractions 23.8%, 18.7%, and 55.4%, respectively). The initial 
TPB density is ~5.2 firrT^. 

126x126x126 computational grid. The governing evolution equations are solved using 
the finite-difference algorithm implemented using Open-MP. No-flux boundary condi- 
tions are imposed on the computational domain boundaries to reflect the assumption 
that materials that comprise the SOFC anode are conserved. 

In the simulated sample, the volume fractions are 23.8%, 18.7%, and 57.5% for 
the Ni, pore, and YSZ phases, respectively. Each phase region is categorized as an 
active, isolated, or dead-end cluster according to the procedure described in Sec. H] By 
examining the YSZ phase with the aforementioned nearest-neighbor scheme, we find 
that, in the initial sample, the entire YSZ phase is fully percolated and active within the 
volume. Using the same procedure, we find that there are 89.1% active, 6.2% isolated, 
and 4.7% dead-end clusters within the Ni phase, and 94.7% active, 1.9% isolated, and 
3.4% dead-end clusters in the pore phase. 

Because our 3D microstructure data are represented as a set of voxels with values 
corresponding to various phases, this categorization procedure seems similar to the 
problem that consists of identifying the percolating clusters in a finite system on the 
basis of a simple cubic network with a coordination number of 6. According to the 
percolation theory, the site percolation of a phase in a simple cubic network is achieved 
when the volume fraction of that phase is above a threshold value of 0.3116. The 
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YSZ phase in our sample is fully percolated and active because its volume fraction is 
much higher than this threshold value. However, it is surprising that nearly 90% of 
the volume of the Ni and pore phases are active while their volume fractions are much 
lower than 0.3 116. This finding may be attributed to two reasons. First, the functional 
layer of our Ni-YSZ anode sample is fabricated by sintering a 50/50wt% NiO-YSZ 
mixture and exposing it to humidified H2 to reduce NiO to Ni. This specific fabrication 
process results in highly percolated pore and Ni phases, even at volume fractions below 
0.3116. Second, in our characterization algorithm, a cluster is considered active if it 
connects any two domain boundaries of the sample volume. Therefore, a cluster with 
a length scale smaller than the sample dimensions can be active as long as the cluster 
connects two neighboring domain boundaries of the sample. In terms of identifying 
the percolated clusters, this criterion is less stringent than the percolation theory that 
requires that the cluster is percolated over an infinite volume. 
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Figure 3: The pore volume change in three different categories over 500 hours of coarsening. The change in 
each category is relative to its initial value. 

A crucial question arises is how these phase clusters evolve with coarsening. The 
coarsening dynamics are greatly simplified in Model B because the Ni surface is the 
only mobile interface in the system and the YSZ phase is immobile. While coarsening 
progresses, the system reduces its total energy by reducing mobile regions with high 
curvatures (which possess smaller length scales and larger surface to volume ratios). 
Because the Ni and pore phases share the same Ni-pore interface, the mobile surfaces 
are the same for both phases. The driving force for material transport thus depends on 
the curvature gradient, which is inversely linked to the characteristic length scale that 
is the inverse of the interfacial area per unit volume of the phase (5^'). For example, 
the characteristic length scale of the Ni structure is defined by the volume of Ni divided 
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by the total interfacial area of Ni. Note that this definition of the characteristic length 
scale is well suited for analyzing multiphase composites with unequal volume fractions 
because it takes into account the effects of volume on the typical size of a phase domain. 
For example, if the volume fraction is large, then we expect that the typical size of the 
domain be larger even if the surface area is the same (this is the case in a two-phase 
system where that phase is the majority phase). Using this definition, an average Ni 
length scale 36% larger than the average pore length scale is found in our sample. 
Inactive pore regions should thus evolve faster and possibly merge with active regions. 

According to the simulation results, during the first 500-hr simulation, the Ni vol- 
ume fractions in each category and the mean Ni length scale remain roughly constant. 
In contrast, as shown in Fig. [3] the dead-end and isolated pore volumes decrease signif- 
icantly, which confirms the faster evolution of regions associated with smaller length 
scales during coarsening. Interestingly, the active-pore volume increases in our sim- 
ulation due to the fact that the smaller isolated and dead-end clusters merge into the 
larger pore clusters, that are most likely in the active cluster category. One interest- 
ing observation is that even though clusters with higher curvatures possess higher free 
energy than those with lower curvatures, they may be trapped in a local equilibrium 
state because the transport pathway is confined to the Ni-pore interface. This kinetic 
constraint explains why the isolated and dead-end clusters in the pore phase did not 
suffer major loss of volume after 500 hours of coarsening. 
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Figure 4: TPB density cliange in each category over 500 hours of coarsening. The active TPB length reduc- 
tion is ~89% of the total TPB reduction after 500 hours. 

The evolution of TPBs is correlated with the evolution of the phases that comprise 
the SOFC anode. The material transport among different clusters of each constituting 
phases not only changes the distribution of the TPBs but also dictates their activity. As 
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shown in Fig. ID the reduction trend of the active TPB density agrees well with the total 
TPB density. The overall TPB density decreased by 20.2% after a 500-hour coarsening, 
while the active and isolated TPB density decreased by 31.4% and 22.3%, respectively. 
In contrast, the unknown TPB density slightly increased by 1.6%. Most of the TPB 
density reduction comes from the 31.4% decrease in the active TPBs. At the early 
stage of coarsening, the rapid reduction of active TPBs is due to the local coarsening 
of high-curvature microstructural features. The small increase in the unknown TPBs 
observed in the early stage of evolution stems from coalescence of inactive Ni domains, 
which occurs rarely. Since the evolution of TPBs is sensitive to local microstructural 
features, it is only indirectly correlated to the coarsening kinetics. 

To summarize, our findings indicate that the coarsening of a Ni-YSZ anode con- 
sumes the mobile phase clusters with high curvatures. The isolated and dead-end pore 
clusters are thus the first clusters that evolve, leading to coarsening of the pore phase, 
which is indicated by the change of 5 j^J^^^^.^^ = 0.1 105 ijm to Sy^^^^^^.^-^ = 0.1 154 yum, i.e., 
a 4.4% increase over 500 hours. While the Ni phase has the same mobile area as the 
pore phase, its mean characteristic length scale, which is defined by Sy^^^^^^y is signifi- 
cantly larger than that of the pore phase and results in a slower evolution. In addition, 
the Ni-pore area reduction is found to be balanced by the increase of the Ni-YSZ inter- 
facial area. The coarsening has thus no significant effect on the mean length scale of 
Ni, which varies from ■SyJ^.j = 0.1711 fim to ^^J^.j = 0.1719 yum (which corresponds 
to a 0.5% change during the course of our simulation. In turn, the coarsening leads 
to a significant reduction in TPBs and the active TPB sites decreases in a very similar 
fashion as the overall TPBs. However, the TPB evolution is sensitive to microstructural 
details and is only weakly linked to the evolution of bulk phases during coarsening. 

5.2. Tortuosity of pores 

Tortuosity represents the geometric aspect of the transport property of a phase that 
comprises a composite. Therefore, the tortuosity of the pore phase plays a crucial role 
in the generation of electricity because the fuels need to be transported through the 
pores to reach the reaction sites. The tortuosity factors of the pore phase are evaluated 
in three orthogonal directions over a 500-hr coarsening period. As shown in Fig. |5] 
although the absolute values differ, the tortuosity factors decrease moderately in all di- 
rections during coarsening. This trend is consistent with the fact that the active pore 
volume increases slightly during coarsening (see Sec. 15.11 Fig.[3]i. These results indi- 
cate a coalescence of pore clusters and a lack of breakup of active clusters, which lead 
to somewhat enhanced transport properties. 

5.3. The active Ni surface 

One advantage of SOFCs is their ability to operate with hydrocarbon fuels as well 
as hydrogen. This is because hydrocarbon fuels can be internally reformed to hydrogen 
with the catalytic effect of Ni at the typical operating temperatures of SOFCs. This 
reforming rate is highly related to the active Ni surface area. Depending on the fuel, the 
reforming kinetics can be very different and involve multiple elementary reaction steps; 
however, the reforming reactions of all types of fuels involve gas species and electron 
transfer, regardless of the detailed reforming mechanisms. Therefore, the active Ni 
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Figure 5: Tortuosity factors in three orthogonal directions of the pore phase over 500 hours of coarsening. 



surface for hydrocarbon reforming must reside at the interface of active pore clusters 
and active Ni clusters. 

To evaluate the catalytic performance of the Ni phase in SOFC anodes, an important 
question to ask is how coarsening affects the active Ni surface area. The evolution of 
an active Ni surface area is dictated by the evolution of active clusters of Ni and pore 
phases. The active Ni surface is identified as the interfaces between active Ni and active 
pore clusters. As shown in Fig. |6] the active Ni surface decrease by 9% after 500-hr 
coarsening, while the total Ni-pore interface decrease by 13.6%. Our finding indicates 
that coarsening reduces more inactive Ni surfaces than those in contact with the active 
clusters of Ni and pore. This is due to the fact that smaller clusters, which have smaller 
length scales and thus have larger driving force toward coarsening, coarsen faster than 
the larger ones. 

5.4. Size effects and evolution kinetics 

We next simulate coarsening in Ni-YSZ anode samples of different sizes, namely, 
(3.2 fim)^, (4.0 fim)^, and (4.8 i-im)^, to examine the size effects. Each of these samples 
is acquired from a portion of the original 914.55 firn^ sample. Specifically, the sample 
of dimensions (3.2 i-im)^ is cropped from the sample of dimensions (4.8 //m)^, while the 
sample of dimensions (4.0 i-im)^ belongs to a different portion of the original sample. 

The Ni surface areas and the TPB lengths are compared over the first 300 hours of 
the coarsening period. The evolution kinetics can be best interpreted from the evolution 
of the Ni surface area in our Model B simulations. As shown in Figs. |7] and |8] the 
Ni surface area reduction rate of the (4.0 fim)^ and the (4.8 fiva)^ specimen are very 
similiar, while the reduction rate of the (3.2 i-im)^ specimen deviates from the other two 
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Figure 6: Ni surface area reduction over 500 hours of coarsening. 
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Figure 7: Reduction of the Ni surface areas for tliree different sample sizes over 300 hours of coarsening. 



significantly even though the specimen is a portion of the (4.8 fim)^ sample. In addition, 
some reduction steps appear in the (3.2 /im)^ curve whereas the other two curves are 
relatively smooth. This suggests that the two larger sample sizes may be sufficient 
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Figure 8: Reduction of the Ni surface area per Ni volume for three different sample sizes over 300 hours of 
coarsening. 



to eliminate the boundary effects from the simulations, while smaller volumes would 
likely suffer from them and may not contain enough particles or statistics for coarsening 
simulations. 

As shown in Fig. |9] the change of the TPB densities of the three specimens are 
compared over 300-hr coarsening. Although the behavior of the TPBs is only indirectly 
related to the evolution kinetics, we find similar trend in the TPB evolution among the 
three specimens as those found in the Ni surface area. Although the TPB density 
reduction kinetics are different, after 300-hour coarsening, the difference of the TPB 
densities between the two larger samples is less than 3%. Since TPB evolution is more 
sensitive to microstructures, this difference may serve as the indication of the local 
variation due to microsctructural details. Thus, in order to determine whether the two 
larger sample volumes are sufficient to represent the microstructure of the anode, we 
would need to simulate a larger sample volume. 

Although rapid reductions of TPBs and the Ni surface area at early stages of coars- 
ening is a result of the large thermodynamic driving force for mass transport, which 
mostly results from regions with high curvatures, the overly steep slopes in Figs. |7] 
and|9]deserve further investigations. 

To examine the length scales involved in the microstructures, the mean curvature 
and the interfacial shape distribution (ISD) of the Ni interfaces in the (4 fJ.m)^ specimen 
are plotted in Fig.[TO](a) and (b), respectively. Due to the diffuse interface nature of 
our models and the model parameters we utilized in our simulations (based on the 
resolution that made simulations feasible), very small microstructural features with 
absolute principal curvature values larger than 0. 17 (which accounts for approximately 
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Figure 9: TPB densities for tliree diflfernt sample sizes over 500 liours of coarsening. The initial TPB 
densities of the (3.2 /im)', (4 //m)', and the (4.8 fim)^ samples are 5.45 fii'rC^, 5.22 fiiiT^, and 5.25 finC^, 
respectively. 

30% of the Ni interfacial area) are not fully resolved. As shown in Fig.[TO](c) and (d), 
if those numerically under-resolved regions are mobile, they will be smoothed out very 
rapidly, leading to coarsening that is faster than expected. 

This numerical smoothing due to under-resolution contributes in part to the rapid 
evolution in the early stage of coarsening simulations and leads to an overestimation of 
the TPB reduction rate. However, these mobile, high-curvature regions inherently pos- 
sess very large thermodynamic driving force and will coarsen sooner or later. There- 
fore, even though the rapid evolution at the early stage of simulations may be caused by 
the insufficient resolution of microstructures, the stabilized value of TPB density after 
further coarsening remains unaffected by it. That is, the prediction of the stabilization 
of TPBs, as well as its predicted value, remains robust, even though the early kinetics 
may be overestimated. Since the microstructures of Ni- YSZ anodes are found to stabi- 
lize after coarsening for a period of time, the final or stablized TPB density rather than 
the short-term evolution kinetics is of importance for the operation of SOFCs. 

6. Conclusion 

We developed two models to study the coarsening kinetics in the Ni-YSZ anode. 
An asymptotic analysis was conducted to Unk our simulation results to the physical 
system. The size eff'ects were studied and an error analysis was performed to validate 
our models. For the model parameters selected, no obvious boundary effects were 
observed when the simulated domain was larger than (4 /im)\ even though the sample 
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size was still insufficient to be statistically representative of the entire microstructure at 
this volume. In addition, the use of the selected mesh resolution had very minor effects 
on the quasi-equilibrium contact angles and evolution kinetics. While the short-term 
evolution kinetics may be affected by the insufficient resolution of microstructures, the 
amount of TPBs and other properties after long-term coarsening can be identified by 
our simulations. 

Although our models contain approximations and the simulation results may be 
affected by uncertainties in the material properties, a reasonable agreement could be 
found for the TPB length reduction between coarsening experiments and simulation 
results using Model B. Unlike Model B, Model A overestimates the TPB reduction due 
to the evolution of the YSZ structure that is induced by an excess mobility at triple 
junctions built into the model. Our simulations show that the major portion of TPB 
reduction occurs during the early stages of coarsening and that stability of TPBs is 
observed provided that YSZ is nearly immobile. 
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The evolution of active TPBs, active Ni surface areas, and tortuosity of the pore 
phase are investigated in our simulations. We found that coarsening has a smaller 
impact on the active Ni surface areas than on their total respective amounts or their 
inactive counterparts. In addition, the tortuosity of the pore phase was found to de- 
crease slowly in all three orthogonal directions during coarsening. These phenomena 
are due to the fact that percolated phase clusters that dictate the active parameters and 
transport have typically larger characteristic length scales, Sy\ than isolated/dead-end 
phase clusters. Smaller clusters experience larger driving forces for coarsening, while 
larger clusters may grow at the expense of these regions if transport can be facilitated 
between them. Active Ni surface areas are also less reduced than their dead/isolated 
counterparts, and the transport property of the pore phase is slightly enhanced. In con- 
trast, the reduction of active TPBs are found to account for most of the total TPBs 
reduction, which suggests that the TPB evolution is sensitive to microstructural details 
and is indirectly related to the evolution of bulk phases during coarsening. 

The proposed coarsening models provide insights into directing the design of anode 
microstructures. The model framework is general and can be applied to many other 
three-phase coarsening systems. Further experiments can help validate and improve 
these models. 



7. Appendix 



7.1. Asymptotic analysis 

The asymptotic analysis that was used is similar to that described in ||43l. Briefly, 
we first determine the outer solution. Substituting equation |22] into |2T] dropping the 
overbars, and replacing the variables 6 with ^, with jj-omC^, t, (), and (p with (pomC^, t, 
we obtain the non-dimensional evolution equations 



dt 



df 
d(pout 



df 



7.1.1. Outer expansion 

An outer expansion is performed by expanding the field i 
in powers of f as follows: 



(27) 
(28) 

and yu„,„ in our model 



cf,„M, t, = t) + cplllcf, tx + 4>2^, t)e + ■ ■ 

fiourC^, t, - t^Zc^, t) + t)^ + pficf, t)e + ■ ■ 

Substituting Eq.|29]into Eq.|27] we obtain to the zeroth order in ( 

( df 



= 0. 



(29) 
(30) 

(31) 
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Eq. [3T|is satisfied with the solution of 

cPZ^^) = VI? e Q_, (32) 

where Q.± are the two bulk phase regions. This solution asserts that far from the inter- 
face, we have one of the equilibrium phases. By substituting Eqs.|29land[3Qlinto Eq.[28l 
and collecting the zeroth order of ^, we find 

= (33) 

Substituting Eq. [32]into Eq. l33]vields 

fiZC^ € fii) = 0. (34) 

7.1.2. Inner expansion 

To denote the inner solution, we replace the variable S with (, /i with /i,„ and (p 
with <pi„. To facilitate the inner expansion, we introduce a moving coordinate system 
(M.C.S.): one coordinate r is parallel to V0, while the other coordinate is the arclength 
s along the interface (which is located at r = 0). The operators close to the interface 
(r ^ 0) are given by: 

V=e,— — , = -^+^,— + -^, (35) 

or OS ar^ or as 

where Kc is the local curvature. Introducing a stretched variable z — rj^, 

V-7-— - ——+K-— + — (36) 

6 dz ds' 6~ dz} 6 dz ds^ 

d I D\ \ d d 

where V„ and are the normal and tangential velocities, and DIDt represents the ma- 
terial derivative, which is the derivative taken based on the M.C.S. Multiplying Eq.lZTl 
by ^2 and rewriting it in terms of a new field variable in the M.C.S. yields 

Expanding the fields in the power of (, gives: 

0,„(z, s, t, = cfP(z, s, t) + 0|^'(z, s, t)^ + cf>f^(z, s, t)e + ... (39) 

f^>.'-fC+l^U + f^i' + - (40) 
Substituting Eq.[39]and|40]into Eq.[38]and collecting the zeroth order of ( leads to: 

,(0)1 



d_ 

d'z 



Mm''' 



dz 



0. (41) 
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Integrating Eq.|4T]with respect to z yields 



a (0) 

dz 



{s,t). 



Taking the limit : 



boo and matching with Eq.[34l we obtain 

a (0) 



(42) 



(43) 



Because (p^^ should vary smoothly from 1 to as z transitions from +oo to -oo, M{(pf^^) 



is expected to be nonzero. We thus have 



94 



(0) 



dz 



= or fif\z,s,t)^gi(s,t). 



(44) 



Matching Eq. |44] with the zeroth order in outer field, the profile of along the z 
direction is governed by 

df d 



= 



dz- ' 

which leads to the following relationship that is useful in change of veriables: 



dz 



Collecting the first order of ^ in Eq. [38] gives: 



d_ 

dz 



dfi 



dz 



,(0) 



dz 



(45) 



(46) 



(47) 



which in turn gives fr.^^' - //^^''(i), or yuj^j' is independent of z- Similarly, collecting the 
second order of ^ in Eq. [38] we find that - l-ifj\s), or fif^ is independent of z- 
Finally, by collecting the third order of ^ in Eq. [38]and considering that yu^^'* and jjf^ 
are independent of z, we have; 



4.(0) 



V„- 



dz 



d_ 

dz 



Integrating Eq. [48]from ; 



dz 

-oo to z = oo gives: 



d 



ds 



/ M{<p'^^)dz 



djJ 



(1) 



ds 



|2„(1) 



d'^l 



ds 



(48) 



(49) 



where we assume a constant I - J M((pf^)dz as it is independent of s. Also, substitut- 

— oo 

ing Eq.[39]and[40]into Eq.[22]and collecting the first order terms of ( leads to: 



(1) 9 f mW 



2^ 52<*(" dd>^^^ 



dz- 



(50) 
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Multiplying Eq.|50]by d(pf^ /dz and integrating over z from -oo to oo gives: 



= -^^ ( I (51) 



Performing a change of variables in Eq.|5T|using the relationship given in Eq.|46lleads 
to: 



1 



^2f(cf>^l^')d<p ^ -K, ■ J, (52) 



I I 

, wherein we assume a constant J - j J2f((pf^)dcf>f^. Combining Eq.|52]and Eq.|49l 



results in: 

V„^I-J-^ (53) 

Performing the integration of / and J (dropping the superscripts and subscripts for 
convenience) leads to: 

1 

7 = J 0(l-0)c/0= (54) 



oo 1 1 _ 

/ = jM((P)dz= J ^^j^d(p^M, J cf>{l-<f>)d<p^^. (55) 

-co ^^■^ 

Substituting I and J into Eq.|53]and restoring the dimensions gives: 
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